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L  INTRODUCTION 


This  manual  documents  the  SMERF  hydrodynamics  code.  SMERF  was  written  at  the  New  Mexico 
Institute  of  Mining  and  Technology,  Socorro,  NM,  and  the  Naval  Air  Warfare  Center  Weapons  Division,  China 
Lake,  CA.  The  code  has  general  applicability,  but  has  been  most  exercised  toward  simulation  of  shock  initiation 
and  sympathetic  detonation  of  explosives. 

SMERF  is  a  combination  of  procedures  used  in  many  well-known  and  established  hydrocodes.  The 
Lagrangian  phase  of  the  calculation  uses  a  modified  Lax-Wendroff  method  (Reference  1).  The  fluxing  of  materials 
in  the  Eulerian  phase  employs  the  van  Leer’s  advection  scheme  (Reference  2)  and  the  Simple  Line  Interface 
Calculation  (SLIC)  algorithm  (Reference  3).  These  are  standard  techniques  that  have  proven  to  be  robust  and 
accurate  numerical  procedures  through  years  of  use  in  the  computational  continuum  dynamics  community.  The 
strength  and  uniqueness  of  SMERF  lies  in  the  detonation  physics  package.  Several  high-explosive  bum  models  are 
available  in  the  code,  and  they  enable  the  user  to  simulate  a  wide  range  of  complex  problems  involving  energetic 
materials.  Sympathetic  detonation  and  shock  initiation  are  examples  of  problems  in  which  SMERF  excels. 

In  Section  II  of  this  manual,  the  differential  equations  of  motion,  along  with  their  finite  difference 
approximations,  are  presented,  and  the  solution  procedure  discussed.  Section  III  describes  inputs  and  outputs.  Two 
example  problems  (shock  initiation  of  high  explosive  by  fragment  impact,  and  sympathetic  detonation  in  a  bomb 
stack)  useful  in  learning  how  to  run  SMERF  are  discussed  in  detail  in  Section  IV. 
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II.  UNDERSTANDING  SMERF 


A.  GOVERNING  EQUATIONS 

The  equations  of  fluid  motion  are  written  in  axis-symmetric  coordinates  below.  An  elastic  constitutive 
model  and  general  functional  forms  for  energetic  reactions  and  equations  of  state  are  also  included. 

Conservation  Equations: 


Mass 

dp  TTdp  dp  fl  drU  dV 

—  +  U  —  +  V-  =  -p - +  — 

dt  dr  dz  r  dr  dz 

Radial  Momentum 

dU  TTdU  dU  ifdiP-SJ  dS„  2S„+S  1 

dt  dr  dz  p\_  dr  dz  r 

Axial  Momentum 

dV  TTdV  rdV  l[d(P-S..)  dS_  SLl 

— +u — +v —  = - - - —+ — -+— 

dt  dr  dz  p[  Sr  dr  r 

Total  Energy 

dE  TTdE  dE  1  fl  drUP  dVP  1 

dt  dr  dz  p\_r  dr  dz 

1  f (  dUSrr  US „  "j  |  (  dVS.,  US7:  ^  |  ( dUS B  |  dVSn  V 
p  v  dr  r  J  v  dz  r  J  {  dz  dr  )_ 

Internal  Energy 

1  =  E  -j(u2  +V2) 


A1 


A2 


A3 


A4 


A5 
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Constitutive  Equations: 


Equation  of  State: 


<®rL  +  U^  +  V^  =  2\x 
dt  dr  dz 


dU  1 
dr  3 


D 


+  2  R'S^ 


^  +  U^-  +  V^-  =  2[i 
dt  dr  dz 


«Z-lD 

dz  3 


-2  R'S, 


dS„  rrdS„  „dS, 


dt 


+  U^  +  V- 


dr  dz 


M- 


dU  dV 
dz  dr 


-R'(S„-S„) 


A6 


A7 


A8 


P  =  P(p,I,W) 


A9 


Energetic  Reaction  Rate: 


dW 

dt 


+  U 


dW 

dr 


+  V 


dW 

dz 


f(P,T,W) 


A10 


Quantities  appearing  in  the  above  equations  are  defined  in  Table  1 . 
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TABLE  1.  Definition  of  Variables  and  Units. 


Variable 

Definition 

Variable 

Definition 

radial  coordinate  (cm) 

P  : 

pressure  (Mb) 

axial  coordinate  (cm) 

deviatoric  stress  (Mb) 

time  coordinate  (ns) 

D 

divergence  (1(ps) 

p 

density  (g/cm3) 

R’ 

spin  tensor  (Kps) 

:  / 

internal  energy  (Mb-cm3/g) 

i§! 

shear  modulus  (Mb) 

E 

total  energy  (Mb-cm3/g) 

w 

burn  fraction  or  porosity 

1 . V 

radial  velocity  (cm/ps) 

T 

temperature  (K) 

■:v 

axial  velocity  (cm/ps) 

The  equations  are  written  in  a  reference  frame  fixed  with  respect  to  the  laboratory  (i.e.,  an  Eulerian  frame). 
In  this  frame,  a  fluid  element  moves  through  space  with  the  local  fluid  velocity  and  changes  because  of  velocity 
gradients.  The  second  term  on  the  left-hand-side  of  the  time-dependent  equations  is  due  to  the  motion  of  the  fluid 
with  respect  to  the  laboratory  frame  and  is  called  the  advection  term.  When  solving  the  equations,  it  is  convenient 
to  treat  these  terms  separately.  First,  the  equations  without  the  advection  term  (Lagrangian  equations)  are  solved. 
The  result  is  what  an  observer  traveling  along  with  a  fluid  element  would  see  and,  therefore,  gives  the  solution  in 
this  moving  frame.  To  get  the  desired  result,  the  Lagrangian  solution  is  mapped  back  to  the  fixed  Eulerian  frame. 
This  mapping  is  called  fluxing  (as  well  as  advection)  because  one  can  easily  picture  it  in  terms  of  fluid  moving 
through  surfaces.  After  the  differential  equations  are  solved,  the  fluid  pressure  can  be  obtained.  It  is  natural, 
therefore,  to  divide  the  solution  procedure  into  three  parts:  1)  Lagrangian  phase  for  source  terms,  2)  Eulerian  phase 
for  the  remap  or  fluxing,  and  3)  pressure  calculation  in  the  equation  of  state. 


B.  NUMERICAL  SOLUTION 

1.  LAGRANGIAN  PHASE 

The  first  step  in  the  solution  procedure  is  to  solve  the  time  dependent  equations  without  the  advection  term. 
This  step  is  called  the  Lagrangian  step,  as  it  provides  the  solution  in  the  moving  Lagrangian  frame.  A  tilde  symbol 
(~)  is  used  to  denote  the  Lagrangian  updated  values. 


Conservation  Equations: 


Bl.l 


PiJ 


s.  pn+'A  _  pn+lA 
„  ot  ri+ytj  ri-y2j 


dr 


B1.2 
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Constitutive  Equations: 


Equation  of  State: 


yn+l  =  yn 

ij  ij 


B1.3 


Fn 

Z'iJ 


St_ 

p  Ij 


{rUP)n+'A 


\»+x 


f+XJ 


3 


"+lA 


r-5  r 


X 


6 z 


B1.4 


[sJl‘ =(sJiJ  +  2»St 


(  tj^Yi  _r/"+>2  i 
I  i+k’J _ i-Y'j  „  _  £)n+Yi 


5  r 


+  2  R'(SK)"j8t  Bi.5 


[sX  =  {S^+2^ 


(  yn+Yi  —  Vn+Y*  i 
I  tJ+Y  iJ~Y  _  _  nw+>2 


5z 


—  A  / 

3  /J 


-2«'(sj;,«/  Bi.6 


rr”+>2  rr»+K  vn*Yi  _  t/w+>2 
+  <+&./  Hk/ 


8z 


8r 


B1.7 


pn+ 1 

'  J 


B1.8 
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Energetic  Reactions: 


B1.9 


Divergence: 


D\ 


,/t+l 

u 


r  II "+^ 
ri+Xui+%J 


■  r  lJn+ % 
r>-Yi  ‘-'AJ 


r;br 


+ 


B1.10 


Material  strength  terms  have  been  omitted  from  the  conservation  equations  above.  These  contributions  are  given 
separately  below. 


un+1  =  un*l+— 

,J  ’-j  p» 

rij 


(Srr)n+v  - {SrrrVv  {snY+Vlv  - {sr)n+y\/  2{srr)n+/l  +  {s„. TYl 

v  rrJl+YiJ  v  v  rzJhj-/2  _J_  v  rrJt,j  v  -JIJ 

6  r  5 z  r; 


Bl.l  1 


y  n+l  _  y„+ 1 

U  U  p" 


v  ^Jij+y2  v  --  .  v  v  ,  v 


6  z 


Sr 


+  - 


B1.12 


En+l=En+ ’+  — 
w  o" 


(KJj"4*  ~{usrr)n+v-  (VSjn+/2v-{VS,X  v  (USrry+/l  +{vs Jn+X2 

_  "J>+YiJ  \  rY_^cl  4.  v  ~  v  "  v  rr/i,y  v  B1  13 

8r  Sz  r, 


+ 


(vsn)n+X2  - {vs^fY2 .  (usj+y2  - {usnY+x 

v  v  2-!zYy  1  ’ELyhA 

5r  5z 
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In  the  differenced  equations,  dr  and  Sz  are  the  horizontal  and  vertical  dimensions  of  a  grid  cell  and  8/  is  the 
time  step.  Subscripts  i  and  j  indicate  the  space  centering  and  superscript  n  the  time  centering.  Quantities  at  cell 
boundaries  and  time  («)  are  weighted  averages.  Quantities  at  the  half-time  step  ( n+1/2 )  and  cell  boundaries  ( i+1/2 ) 
and  (j±l/2)  are  computed  as  follows: 


yn+Yi  -y» 
yiJ+Yz  y'J+ 


Pn+Yi  -  pn  _^foc2)" 

'+><,/  “  r‘*YiJ  2  » 


'i+YiJ 


bt  12 

~P"  _ 

r/+lJ 

P'+X.r 

5r 

8t/2 

’  pb  _ 

<j+i 

P "j.K 

5z 

•wu 

l,J 


-Yi 


r*y?r 


6z 


B1.14 


B1.15 


B1.16 


pn+Yi  -  pn  _^1(DC2)" 

u+Yi  u+Yi  2  P  u- 


<J+Yt 


Mae* 


YuX-yrJ,yt  ,  r^-v-j 

rfir  5 z 


B1.17 


The  half-step  pressure  update  has  used  dP/dt=-pc^D,  a  prognostic  equation  which  can  be  derived  using  Equations 
A1-A3.  The  half-step  stresses  are  computed  as  follows. 


(S,Xl=(S,XKj-2rft 


Uw-U’u  If 

5 r  3 


KrvY i  iD. 

Sr  3  ■  K 


(SXA-(s^h,y,-^1 


trm-K  lD. 


5  r 


5 z 


1+1/2, 7 


Ir. 

Sr  +  Sz  3 


B1.18 


B1.19 


B1.20 


B1.21 
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2.  EULERIAN  PHASE 

The  second  step  in  the  solution  procedure  is  to  map  the  solution  obtained  in  step  1  (Lagrangian  update) 
back  onto  the  fixed  Eulerian  mesh.  This  is  accomplished  by  moving  material  across  zone  boundaries.  A  one¬ 
dimensional  algorithm  is  used  for  the  flux  calculations,  with  each  direction  treated  separately.  Consider  the 
movement  of  mass.  Let  5x  be  the  flux  volume  swept  out  by  moving  the  grid  line  at  the  right  edge  of  cell  (idj)  back 
to  the  fixed  Eulerian  position.  We  use  an  upstream  scheme  in  which  the  donor  cell  value  of  p  is  evaluated  at  the 
center  of  the  flux  volume  using  one  of  three  candidate  slopes  (*S7,  Sc,  Sr).  These  concepts  are  illustrated  in  Figure  1 . 


Eulerian  Position 


s/\;- 

\  Lagrangian  Position 

/ 

"4  if 

Sc  '  •  '  . 

1 

'\Sr 

Flux 

Volume 

1 

id-1 

id 

id+1 

spatial  coordinate  (r) 


FIGURE  1.  Illustration  of  van  Leer's  Advection  Scheme. 


The  flux  volume  is 


St 


B2.1 


The  mass  in  5x  crossing  the  i+1/2  boundary  is 


8m, 


i+YiJ 


B2.2 


where  the  density  p  is  evaluated  at  the  center  of  the  flux  volume  using  the  smallest  of  three  candidate  slopes 


g  _  P idj  Pid-\J  g  Pid+\J  Pid-lJ  g  Pid+\J  P idj 

l~  5r/2  ’  c~  25 r  ’  r  _  8r/2 


B2.3 
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dp\  MAT(|5j,|5e|,|5,|)*5/GA^(S'/  *Sr)  if  S,*Sr>  0 
dr  |0  if  S,*Sr<  0 


The  new  (Eulerian  updated)  mass  in  the  cell  is  simply  the  net  flux, 


<  =  +&%*/■ 


'  dm, 


‘-YiJ 


B2.5 


The  fluxing  of  other  quantities  is  handled  in  the  same  way  as  the  mass.  Letting  XP  represent  the  momentum 
or  total  energy  and  vj/  the  mass  specific  values  (velocity  or  specific  total  energy),  the  amount  of  lP  crossing  a  cell 
boundary  is  n.  Where  vy,  like  p,  is  evaluated  at  the  center  of  the  flux  volume.  Conservation  of  ¥  then  gives 

the  new  (Eulerian  updated)  values  of  \\i. 


<'  = 


B2.6 


The  fluxing  scheme  described  above  is  van  Leer’s  modification  (Reference  2)  of  the  upstream  advection 
scheme.  The  algorithm  is  second-order  accurate  and  does  not  introduce  spurious  oscillations  common  to  higher 
order  schemes.  When  more  than  one  material  is  present  in  the  donor  cell,  its  kind  and  the  order  in  which  the  flux 
volume  is  filled  must  be  determined.  The  SLIC  (Simple-Line-Interface-Calculation)  algorithm  (Reference  3)  is 
used  for  this  purpose.  SLIC  looks  in  the  zone  ahead  (downstream)  and  behind  (upstream)  the  donor  cell  and  assigns 
one  of  four  index  number  pairs  to  each  material  in  the  donor  cell.  The  pair  (0,0)  means  that  none  of  this  particular 
material  is  present  in  the  zone  behind  the  donor  cell  and  none  is  present  in  the  zone  ahead  of  the  donor  cell. 
Accordingly,  (0,1)  means  that  there  is  none  behind  and  some  ahead,  (1,0)  indicates  some  behind  and  none  ahead, 
and  (1,1)  says  that  there  is  some  of  this  material  in  the  zones  both  upstream  and  downstream  of  the  donor  cell.  For 
positive  velocities,  the  materials  are  fluxed  in  the  order  [(0,1), (1,1), (0,0), (1,0)]  and,  for  negative  velocities,  the  order 
is  [(1,0), (1,1), (0,0), (0,1)].  Figure  2  depicts  the  algorithm  for  a  simple  case  in  which  material  flows  from  the  center 
cell  containing  two  different  materials  into  the  cell  on  the  right.  A  two-dimensional  example  of  how  SLIC 
represents  a  fluid  interface  is  shown  in  Figure  3. 


actual  fluid  interface  SLIC  representation  SLIC  representation 

horizontal  pass  vertical  pass 

FIGURE  2.  Illustration  of  How  the  SLIC  Algorithm  Represents  an  Interface  for 
Movement  in  Two  Directions. 
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FIGURE  3.  Simple  Illustration  of  How  the  SLIC  Algorithm  Assigns  an  Order  to  the 
Movement  of  Material  Across  a  Cell  Boundary. 


C.  EQUATIONS  OF  STATE 

L  ANALYTIC  FORMS 

Step  3  in  the  solution  procedure  consist  in  calculating  the  pressure  from  the  newly  determined  density  and 
internal  energy.  Four  equations  of  state  are  available:  Perfect  Gas,  Gruneisen,  HOM  (Reference  4),  and  JWL 
(Reference  5).  In  the  following  equations,  Hugoniot  and  isentrope  reference  curves  are  indicated  by  subscripts  “h  ” 
and  “i,  ”  respectively.  The  compression  is  |i=p/p0-l,  and  e0  is  the  heat  of  reaction  normalized  by  the  initial  specific 
volume. 


Perfect  Gas: 


Gruneisen: 


HOM  Solid: 


Cl. 2 


P  = 


P„C2n[l  +  (l-y/2)n- a\j}  / 2] 


+  (y  +a\i)I  \x>\  C1.3a 


P=p„cli  +  (l,+‘v)1 
P=Pt+jF(i-iJ  rjr>  1 


fO.  <  1  Cl. 3b 


Cl. 4a 
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.  c’fc-r) 

‘  [K-s(K-vjf 

T=Tt+(l-It)/C, 

In  Th  =  F  +  a(lnV)+  H(biVf  +  l(lnVf  +  J(\nVf 


Cl. 4b 


Cl. 4c 


C1.4d 

C1.4e 

C1.4f 


HOM  Gas: 


JWL: 


P  =  P,+±(/-I,)/V  C1.5a 

T=T,+(l-I,)/C¥  Cl. 5b 

ln/>  =  A  +  B{\nV)  +  C{\nV)2  +  Z>(lnF)3  +E{  InF)4  C1.5c 

In?;  =  Q+  R(lnV)  +  S(\nV)2  +  TflnF)3  +£/(lnF)4  Cl.5d 

ln(/,  +Z)  =  K  +  l(ln  P)  +  M(  In  P)2  +  tf(ln  Pf  +  0(ln  Pf  Cl  .5e 

-l/p  =  i?  +  25(lnF)  +  3r(lnF)2+4[/(lnF)3  C1.5f 


P=A 


1- 


vVo 

r,v  j 


-R,VIV„ 


+  B 


1- 


wiy 

R2Vj 


,-R2V/V° 


+  w\ 


I 

—  +  eo 

V% 


T  = 


A  —RV  $  -R  V  I 

—  e  1  + —  e  2  +  —  +  e„ 


V*. 


R, 


v„ 


jC, 


Cl. 6a 


Cl. 6b 


2.  ADIABATIC  PRESSURE  EQUILIBRIUM  IN  MIXED  CELLS 

Some  procedure  for  calculating  the  pressure  in  computational  cells  containing  different  kinds  of  materials 
is  needed  because  the  partial  pressures  are  generally  not  equal  after  the  Lagrangian  and  Eulerian  steps.  We  use  a 
Newton-Raphson  iterative  procedure,  adjusting  the  partial  volumes  (t,)  along  adiabats,  until  the  partial  pressures 
(Pf)  are  equal.  Consider  the  simple  case  of  two  materials.  The  mass  (M=M/+A/2),  volume  (t=t1+t2),  and  internal 
energy  (I=MIII+M2I2)  of  the  cell  as  a  whole  are  known  and  serve  as  constraints  during  the  iteration.  Requiring  the 
partial  pressures  equal  to  the  first  order  on  the  k+1  iteration  step  gives 
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pk+l  _  pk  ^1  _  pk+ 1  _  pk  _j_  dl\  dXg 

1  ~  1  dVx  Mx~  2  2  dV2  M2 

where  the  derivatives  are  taken  along  adiabats 

dP  _dP  SP  dl  _  8P  _pdP 

dV  ~  dV  +  dl  dV  ~  dV  dl 


Using  the  volume  constraint  ( dx ?  =  —  dxx ),  we  solve  for  the  change  in  partial  volumes. 


dx1=(P2k-P1k/ 


J _dPL 
Mx  dVx 


1  dP2 ' 
M2  dV2  , 


CIA 


Cl  2 


C2.3 


The  accompanying  energy  changes  along  the  adiabats  are  and  M2<i/2=-/?2^T2.  Since  this  procedure 

does  not  conserve  the  total  energy  in  the  cell,  the  cell  pressure  P*  l=pk+l=pk+l  defined  by  C2.1  is  used  as  follows 


Jk+ 1  _  Jk  _  pk+¥t 
'  1  M\ 

CIA 

Jk+ 1  —  Jk  —  pk+) 2 

2  2  M, 

C2.5 

which  conserves  the  energy  during  each  step  of  the  iteration.  The  pressure  computation  is  done  after  the  Lagrangian 
step  and  after  the  Eulerian  step.  On  entering  the  mixed  cell  routine,  the  sum  of  the  initial  material  volumes  and 
energies  does  not  equal  the  cell  total  volume  and  energy.  Therefore,  prior  to  the  pressure  equilibration  calculation, 
the  excess  volume  is  distributed  to  the  materials  on  a  volume  basis,  and  the  excess  energy  is  distributes  on  a  mass 
basis. 


3,  THERMODYNAMIC  EQUILIBRIUM  FOR  EXPLOSIVES 

Some  procedure  is  required  for  calculating  the  pressure  in  computational  cells  containing  a  mixture  of 
reactants  and  products.  The  method  we  use  is  to  bring  the  constituents  to  pressure  and  temperature  equilibrium  by 
adjusting  the  partial  volumes  and  energies,  while  conserving  the  total  volume  and  energy  in  the  cell.  Let  us 
consider  the  case  where  the  mixed  cell  contains  unreacted  solid  (s)  and  reaction  product  gas  (g).  We  have  the  solid 
and  gas  equations  of  state 


f,  =  P,(K,I.)  C3.1 

e,  =  r,Ki,)  C3-2 

The  specific  volume  and  internal  energy  are  partitioned  into  unreacted  solid  and  gaseous  products  by  the 
conservation  equations  and  the  unbumed  fraction  W. 


V  =  WVs  +  {\-W)Vg  C3.3 

I  =  WIs  +  (\-W)Ig  C3.4 

We  assume  the  solid  and  gas  phases  are  in  pressure  and  temperature  equilibrium. 
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P  =  P  C3.5 

*  *g 

T  =  r  C3.6 

*  g 


Introduction  of  the  temperature  necessitates  a  thermal  equation  of  state,  which  we  take  to  be  linear  in  the 
temperature 


i,  =  L(K)+c„(vX 

C3.7 

U.* 

if 

+ 

II 

C3.8 

Substituting  C3.7  and  C3.8  into  C3.4  and  imposing  C3.6  yields 

I  =  h+CVT 

C3.9 

where 

i.-wiSr,)+Q-w)iJyt) 

C3.10 

c,=wcjv,)+(i-<ncv{rt) 

C3.ll 

An  equation  for  the  temperature  is  obtained  from  C3.9  and  C3.3 

t  =  t(v,i,w,vs) 

C3.12 

Every  quantity  is  known,  except  for  Vg.  Once  the  equilibrium  temperature  is  known,  equations  C3.7  and  C3.8  give 
the  partial  energies,  Ig  and  IS7  and  then,  finally,  C3.1  and  C3.2  give  the  partial  pressures  Ps  and  P s  One  gets  a 
relation  for  the  pressure  difference  dP  as  a  function  of  the  unknown  gas  volume  Vg  and  there  is  a  particular  choice  of 
Vg  for  which  the  pressure  equilibrium  ( dP=0 )  holds.  Newton’s  method  is  used  to  solve  for  V 


D.  HIGH  EXPLOSIVE  BURN  MODELS 

An  equation  describing  the  rate  of  decomposition  of  energetic  materials  is  required  for  the  simulation  of 
detonation  phenomena.  Four  bum  models  are  available:  Programmed,  Arrhenius,  Nucleation  &  Growth,  and  Forest 
Fire. 


L  PROGRAMMED  BURN 

This  model  is  used  to  generate  a  detonation  wave  in  an  explosive  which  travels  at  the  Chapman-Jouget 
velocity  (DCJ).  The  detonation  is  started  at  a  particular  location  which  is  specified  in  the  input.  The  detonation  is 
assumed  to  start  at  time  zero.  At  the  proper  time,  the  explosive  reactants  within  a  cell  are  converted  to  products. 


2.  ARRHENIUS  BURN 

The  Arrhenius  rate 
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dW 

dt 


involves  two  parameters:  the  activation  energy  (£*),  and  a  frequency  factor  (Z*). 


D2.1 


5.  NUCLEA TION  &  GROWTH 

The  Nucleation  and  Growth  bum  has  the  most  complex  form  of  all  bum  models  in  the  code.  The  bum  rate 
consists  of  three  terms:  an  ignition  term  and  two  growth  terms. 


dW 

dt 


dW 

dt 


dW) 


dt 


D3.1 


Each  term  can  be  switched  on  or  off  depending  on  the  product  mass  fraction.  The  ignition  rate  requires 
that  the  compression  of  the  material  be  greater  than  a  threshold  value  and  that  the  rate  vanish  after  some  fraction  of 
the  material  is  burned.  The  individual  rates  are 


dW\  _ 

oc[(K/v-i')-ccr] 

*k 

0 

dw)  -< 

G,W'a'{l-Wf  Pm 

0 

dw)  -\ 

\G2W^{\-Wf  Pn 

*)*  \ 

1° 

(1-(F)>F_ 

(!-»')<■?«, 

(1  -W)>F^ 

(\-V)<Fmr 

(1  -W)>Fmgr 

(1  -<r)<F.,!r 


D3.2 


D3.3 


D3.4 


The  14  parameters  in  the  Nucleation  and  Growth  model  are  determined  by  comparing  experimental  data 
with  theoretical  calculations.  The  best  experiments  for  this  purpose  use  embedded  gauges  which  record  the  response 
of  the  explosive  to  planar  shocks.  Model  parameters  are  then  adjusted  to  obtain  agreement  between  gauge  records 
and  calculations. 


The  Nucleation  and  Growth  bum  model  originated  with  Lee  and  Tarver  (Reference  6)  and  was  expanded 
by  Tarver  and  Hallquist  (Reference  7)  to  its  present  form.  In  general,  the  values  for  the  nucleation  and  growth 
constants  depend  on  the  choice  of  equation  of  state  for  the  reactants  and  products  and  the  rule  governing  their 
mixtures.  In  previous  work,  the  JWL  equation  of  state  was  used  for  both  phases.  However,  the  mixture  rule  was 
used,  not  the  equilibrium  temperature  and  pressure  rule  used  in  SMERF.  Therefore,  nucleation  and  growth 
constants  found  in  the  literature  cannot  be  used  in  SMERF. 
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4.  FOREST-FIRE 

The  Forest  Fire  bum  model  is  a  global  reaction  kinetics  model  for  the  decomposition  of  heterogeneous 
explosives.  The  bum  model  is  calibrated  directly  using  the  results  of  wedge  tests.  Analytic  methods  for  calibrating 
the  model  were  derived  by  Forest  (Reference  8),  the  originator  of  the  model.  Lundstrom  (Reference  9)  has 
reviewed  and  extended  Forest’s  original  work. 

In  a  wedge  test,  a  shock  wave  is  introduced  into  the  test  explosive.  Due  to  burning  in  the  shocked 
explosive,  the  shock  wave  accelerates  and  transitions  into  a  detonation  wave.  The  progress  of  the  accelerating 
shock  wave  is  monitored  experimentally  using  the  wedge  shape  of  the  explosive.  The  distance  that  the  shock  wave 
takes  to  transition  to  detonation  is  called  the  run  distance.  The  experimental  data  is  customarily  displayed  as  a  log- 
log  plot  of  the  run  distance  (. X *)  as  a  function  of  the  initial  shock  pressure  ( P ).  This  representation  of  the  data  is 
known  as  a  Pop  Plot,  after  its  originator  A.  Popolato  (Reference  10).  For  the  vast  majority  of  explosives,  data  on 
the  Pop  Plot  can  be  fit  very  well  with  a  straight  line,  in  which  case  the  wedge  test  results  can  be  expressed  in  the 
form 


where  S  is  the  slope  of  the  line  and  is  the  pressure  (Kb)  required  to  give  a  one-centimeter  run  distance.  The  result 
of  each  wedge  test  is  a  shock  trajectory  from  which  the  shock  pressure  as  a  function  of  the  shock  propagation 
distance  and  time  can  be  obtained.  An  important  assumption  in  Forest's  derivation  is  the  "single  curve  buildup 
principle".  According  to  this  principle,  all  of  the  wedge  test  shock  trajectories  corresponding  to  different  initial 
shock  pressures  can  be  reduced  to  a  single  curve.  The  trajectories  must  first  be  shifted  in  time  and  space  to 
superimpose  their  points  of  transition  to  detonation.  Equation  D4.1,  therefore,  describes  the  resulting  single 
trajectory. 

In  Forest’s  method,  as  described  by  Lundstrom  (Reference  9),  one  calculates  the  bum  rate  behind  the  shock 
wave  which  is  accelerating  according  to  equation  D4. 1 .  To  do  this,  one  uses  the  differential  equations  expressing 
the  conservation  of  mass,  momentum,  and  energy  in  the  reacting  flow  field  and  the  Hugoniot  equations  expressing 
the  same  conservation  principles  across  the  shock  wave.  An  expression  can  be  derived  for  the  bum  rate  behind  the 
shock  wave  having  the  functional  form 

^  =  F(r,W)-^  +  G(P,W)f-  D4.2 

dt  dX  dx 

where  dP/dX*  is  the  change  of  shock  pressure  with  run  distance  as  given  by  equation  D4.1,  and  partial  derivative  is 
the  pressure  gradient  immediately  behind  the  shock  front.  The  expressions  F  and  G  are  known  functions  of  the 
shock  pressure  (P)  and  the  reactant  mass  fraction  ( W)  behind  the  shock.  In  general,  one  can  assume  explosive 
decomposition  within  the  shock  wave  so  that  W  can  be  less  than  one.  One  normally  makes  the  simplifying 
assumption  that  dP/dX-0,  since  no  information  regarding  this  quantity  is  obtained  from  the  wedge  test.  Together 
with  equation  D4.1  for  dP/dX *,  equation  D4.2  reduces  to 

—  =  W(P,W)  D4.3 

dt  K 

That  is,  the  bum  rate  required  to  accelerate  the  shock  wave  is  a  known  function  of  shock  pressure  and  degree  of 
reaction  within  the  shock  wave.  This  function  can  then  be  used  to  calibrate  different  bum  models.  Forest  assumed 
that  the  reaction  of  the  explosive  takes  place  according  to  a  first-order  decomposition  reaction.  First-order  reactions 
rates  have  the  form 
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1  dW 
W  dt 


=  /(i5) 


D4.4 


where  the  reaction  rate  function,/^,  is  independent  of  the  solid  mass  fraction  and  is  expressed  solely  as  a  function 
of  the  local  pressure. 


An  alternative  to  the  first-order  hypothesis  is  the  ignition  and  growth  concept  used  by  Lee  and  Tarver  to 
describe  explosive  bum.  The  shock  initiation  of  the  heterogeneous  solid  explosive  is  modeled  by  ignition  at 
localized  hot  spots  followed  by  hole  burning  at  the  growing  hot  spot  boundaries.  This  model  has  had  some  success 
in  correlating  the  detailed  mass,  velocity,  and  pressure  histories  behind  the  shock  wave  during  shock-to-detonation 
experiments,  lending  credibility  to  the  basic  concept.  According  to  the  ignition  and  growth  model,  the  global 
reaction  rate  is  a  minimum  at  the  shock  boundary  where  W=1  because  of  the  small  bum  surface  area  of  the  hot 
spots.  This  is  in  direct  opposition  to  the  assumption  of  first-order  kinetics,  where  the  reaction  is  maximum  at  the 
shock  front. 

Lundstrom  (Reference  9)  presented  an  alternative  zero-order  kinetics  assumption  where  the  explosive 
decomposition  is  given  by  the  form 


dW 

dt 


=ap) 


D4.5 


which  ignores  ail  dependence  of  the  reaction  rate  on  the  solid  mass  fraction.  This  model  should  lie  somewhere 
between  the  extremes  of  the  ignition  and  growth  model  and  the  first-order  kinetics  employed  by  Forest.  Lundstrom 
discussed  the  calibration  for  both  first  and  zero-order  reaction  rate  kinetics  using  different  assumptions  on  the 
degree  of  reaction  within  the  shock  wave.  The  resulting  bum  rate  models  were  tested  by  performing  numerical 
simulations  of  the  wedge  tests  to  see  if  the  overall  wedge  test  results,  equation  D4.1,  are  obtained.  It  was  found  that 
the  zero-order  reaction  rate  law  gave  reasonable  results  when  calibrated  using  the  assumption  of  no  reaction  with  the 
shock  wave.  In  particular,  it  was  found  that  the  reaction  rate  was  well  behaved  at  pressures  beyond  the  Von 
Neumann  spike  pressure.  The  original  Forest  Fire  model  became  infinitely  near  the  Chapman-Jouget  pressure.  This 
is  important  when  trying  to  apply  the  bum  model  to  problems  involving  detonation  failure.  For  an  inert  shock,  the 
explicit  form  of  equation  D4.3  is 


where 


dW 

dt 


/(p)=-jl+ 

f  PP\  ~  Py  y 

v  PP\  ~2PV  J- 

n  dP 

„  dP 

P  =  — 

Py  =  - 

1  dl 

r  dV 

Us  dP 

tj? 


4= 


dP 

dW 


D4.6 


D4.7 


and  Us  is  the  shock  velocity.  Given  an  equation  of  state  for  the  explosive  reactants  and  products,  the  Hugoniot 
relations  for  the  shock,  and  equation  D4.1  for  the  run  distance,  the  right  hand  side  of  D4.6  can  be  evaluated  as  a 
known  function  of  the  shock  pressure. 


The  SMERF  code  uses  the  zero-order  kinetics  version  of  Forest  Fire  calibrated  using  the  inert  shock 
assumption.  The  function  f(P)  is  evaluated  at  a  number  of  points  at  the  start  of  the  SMERF  calculation.  They  are 
used  as  a  basis  for  interpolation  using  a  cubic  spline  fit  of  log(j)  as  a  function  of  log(P). 
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E.  POROUS  MATERIALS 

The  porous  material  model  in  SMERF  is  based  on  the  P-a  model  originally  proposed  by  Herrmann 
(Reference  1 1).  A  similar  model  is  incorporated  into  the  CTH  code  in  a  substantially  different  form  (Reference  12). 
A  distention  parameter,  a,  relates  the  overall  specific  volume  of  the  porous  material,  V,  to  the  specific  volume  of  the 
solid  matrix,  Vs. 


V 


Then  the  equation  of  state  of  the  porous  material  is  related  to  the  solid  matrix  material  by 

P(V,I,a)  =  -P(Vs,Is) 
a 


El 


E2 


T(V,I,a)  =  T(Vs,Ii.)  E3 

A  diagram  of  the  equation  of  state  of  the  porous  material  on  the  P-V  plane  is  shown  in  Figure  4.  The 
curves  are  shown  for  constant  internal  energy.  The  curve  labeled  <x=l  represents  the  compression  curve  for  the 
solid  matrix  material.  A  curve  labeled  for  constant  a>l  represents  solutions  of  E2  for  the  distended  porous  material. 


FIGURE  4.  The  Equation  of  State  of  a  Porous  Material  in  the  P-V 
Plane  With  Constant  Internal  Energy. 
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Compression  of  the  porous  material  is  elastic  and  reversible  at  low  pressure.  The  elastic  compression  curves 
can  be  described  by  equation  El  with  a  given  as  a  function  of  pressure  by 


a=f-g(PJ)  E4 

where 

g(0,/)  =  l  E5 

and 

g(P,  1)  =  1  E6 


According  to  E3  and  E4,  the  parameter /  is,  therefore,  the  value  of  a  at  zero  pressure.  It  is  the  extent  of 
irreversible  crushing  that  the  porous  material  has  undergone.  The  parameter,  f  is  a  constant  during  the  elastic 
compression  of  material.  The  shape  of  the  elastic  compression  curve  is  given  by  the  function  g(P,f)  which  will  be 
specified  later.  Several  elastic  compression  curves  of  constant  /  are  included  in  Figure  4.  The  condition  on  the 
function  g ,  Equation  E5,  ensures  that  the  compression  of  the  solid  matrix  material,  a=l,  is  a  special  case  of  E3. 
This  factoring  of  a  into  / and  g  is  not  the  usual  procedure,  but  it  makes  the  geometry  of  the  porous  equation  of  state 
more  physically  intuitive  (at  least  to  the  authors). 

The  porous  material  undergoes  reversible,  elastic  compression  until  it  reaches  a  limiting  pressure  designated 
by  PJKI),  where  the  compression  becomes  irreversible  crushing.  This  crushing  curve  is  included  in  Figure  4.  The 
compression,  therefore,  first  follows  a  constant  /  line  until  it  intersects  the  Pm  line,  which  it  then  follows  until  it 
intersects  the  solid  matrix  line  a=l.  Unloading  follows  the  constant  /  line  which  intersects  the  Pm  line  at  the 
maximum  compression. 

Experimentally,  the  solid  matrix  equation  of  state  is  often  available  or  can  be  measured.  The  irreversible 
compaction  curve,  Pm(V,I),  can  be  measured  in  quasi-static  compaction  tests,  and  the  data  is  frequently  available  in 
a  tabular  or  a  fitted  functional  form.  Shock  testing  of  the  material  can  be  performed  to  give  the  elastic  shock 
precursor  velocity  that  can  be  used  to  determine  the  g(P,f)  function.  The  precursor  velocity  is  related  to  the  slope  of 
the  constant / curve  at  zero  pressure  in  Figure  4. 

The  constant /lines  satisfy  the  implicit  relation  which  results  from  combining  El,  E2,  and  E4, 


f-g{pj)-p  = 


v 


E  6 


Partial  derivatives  of  the  pressure  P  with  respect  to  V  and  E  are  readily  calculated  as 


/, 


E7 


and 
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dP\  dP. 


at,—- 

dl)v  dl, 


E8 


sy  v 


where  the  quantity  is  given  by 


5=i+ 


v  2  dPs 
p+— — - 
V  dV. 


/,  J 


a  lng 


The  sound  speed  satisfies  the  thermodynamic  identity 


a 2  =  -V‘ 


dP 


dP 


dP 


E9 


P— 


[  dVJ!  dlJy  J 

so  that,  for  the  porous  material,  the  sound  speed  is  given  by 

2  1  2 
a-  =  — - a j 

£,a? 

where  as  is  the  sound  speed  in  the  solid  matrix  material.  The  specific  form  for  the  function  g  is  chosen  to  be 

2(f  P)  = - 1 - 

1  +  hP{f  -1) 

where  h  is  a  constant  that  is  chosen  to  match  the  measured  elastic  sound  speed,  ap  at  porosity^,. 

At  P=0,  then  g=l,  a=l ,  and  a=  -  V2S  dPs/dVs,  so  that  Equation  E9  for  £  becomes 


E10 


Ell 


E12 


l(P  =  0)  =  \  +  ha*(f-\)/V 
Substituting  E13  for  £,  in  El  1)  and  solving  for  h,  one  gets 


h  =  - 


V 


(/-l) 


J _ 1_ 

a 2  al 


E13 


E14 


The  calibration  condition,  sound  speed  a=ap  at  porosity,  fp  determines  h.  The  porous  model  is  treated  in 
SMERF  as  an  equation  of  state  with  an  internal  degree  of  freedom,/  A  subroutine  is  required  that  evaluates 

P=P(V,I,f)  and  T=  T(V,I,f)  E15 

To  do  this  in  principle,  P  is  evaluated  as  if  the  material  were  on  the  elastic,  reversible  curve  corresponding  to  the 
input  compaction/  If  the  resulting  pressure  is  greater  than  the  crushing  pressure,  Pm,  then  a  new  compaction, /=/m, 
is  computed  which  satisfies 


Pm(V,I)=P(V,I,fm) 


E16 
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The  subroutine  then  returns  P,  T,  and  the  compaction,/  which  may  be  modified  to  fm  if  required.  In  practice, 
one  first  finds  out  which  region  the  material  is  in  (elastic,  crushing,  or  fully  crushed)  and  computes  the  crush 
pressure,  Pm  =  Pm(V7I).  One  can  find  out  whether  the  material  is  fully  crushed  if  the  following  inequality  holds 

Pm<Ps{V,I)  E17 

One  sets /=/  and  uses  the  equation  of  state  for  the  fully  crushed  material.  If  not,  then  one  estimates  g  as  gm  = 
g(Pm,fj  and  computes  a  trial  pressure  P„ 


P,  = 


1 


'  fs, r '• 


1  v_^ 

yfSm  ’  7 


If 


then  the  material  is  crushing  and  one  must  solve 


P  <P 

1  m  ^  J  l 


P  = 


1 


fmg(Pm,f) 


V 


JnSiPm’/y1  J 


E18 


E19 


E20 


for  fm.  A  non-linear  second-order  solver  is  used  for  this  purpose.  The  solution,  fm  for  the  crushing  material  is 
known  to  be  bounded  by 


/>/->! 

Finally,  if  the  inequality  E19  is  not  satisfied,  then  the  material  lies  on  the  elastic  curve  and  one  must  solve 


E21 


— l- — p 

fg(P,f)  s 


v 


.fg(PJ) 

for  the  pressure  P.  The  same  non-linear  solver  is  used  in  this  step. 


j 


E22 


F.  MATERIAL  STRENGTH 


An  elastic,  perfectly  plastic,  constitutive  model  is  used  to  describe  strength  effects  in  solids.  For  linear 
elasticity  (Reference  13)  the  rate  of  change  of  the  deviatoric  stresses  is  proportional  to  the  traceless  rate  of  strain. 


4p  =  2psap  =  2(x(sap-8rySap) 

Here,  p  is  the  shear  modulus,  and  the  strain  rates  are  defined  by 


Sa|3  - 


rdU„  dUt ' 
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dx  „  dx. 
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FI 
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The  time  derivative  of  the  stress  in  Equation  FI  is  not  objective;  that  is,  the  material  response  will  change  if 
the  underlying  frame  with  respect  to  which  the  stresses  and  deformations  are  measured  undergoes  a  rigid  translation 
or  rotation.  A  variety  of  objective  stresses  have  been  formulated.  The  Jaumann  rate  is  the  most  widely  used  in 
codes.  Hermann  (Reference  14)  examines  the  relative  merits  of  several  stress  rates.  With  the  Jaumann  rate,  our 
constitutive  relation  is 


^ap  ^ay^|3y  ^yp-^xy  “  ^M^aP  “  ^fx(sap  £yySap) 

where  R  is  the  rotation  rate 


F3 
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^ocp  " 


dU  dUn 


dxn  dx, 
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The  individual  components  of  F3  are: 


dt 

dS„ 

dt 


=  2jLl 


feu 

v  Sr 


+  2  R'S^ 
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dSn  ( 8U 

— —  =  u  — 

dt  \  dz 


and  in  two-dimensional  axis-symmetric  coordinates, 


dV 

dr 
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r  dr  dz 

The  plastic  flow  regime  is  determined  by  the  von-Mieses  criterion  when  the  second  stress  invariant,  J2=Sa p  Sap, 
exceeds  the  known  flow  stress,  Y0.  The  individual  deviators  are  then  brought  back  to  the  flow  surface 


Jap 


F10 
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III.  RUNNING  SMERF 


A.  BASIC  INPUT  (The  setup  File) 

There  are  six  input  files  which  SMERF  requires  for  problem  initialization:  setup,  mshgen,  matcor, 
outdata,  gaugein,  and  an  equation  of  state  data  file.  The  entries  in  all  of  the  input  files  have  the  FORTRAN-free 
format  where  the  data  are  separated  by  a  comma  or  a  space.  The  input  files  are  prepared  using  any  convenient  text 
editor.  Basic  problem  options  are  specified  in  file  setup.  One  example  problem  discussed  in  Section  IV  is 
sympathetic  detonation  in  a  stack  of  four  bombs.  For  this  problem,  the  setup  file  looks  like  the  one  shown  in  Table 
2. 


The  variable  names  appearing  within  the  square  brackets  in  Table  2  are  for  convenience  and  a  memory  aid 
and  are  ignored  by  the  SMERF  program.  Entries  in  the  setup  file  are  described  line-by-line  below. 


TABLE  2.  Setup  Input  File  for  the  Bomb  Stack  Example  Problem. 


'Bomb  Stack' 

0,  'pi',  0.0,  0.25 
F.'restartfile' 

1.0, 1.0 

200.0,  1200,  0.4 
F,  .2,  .3,  100,0, 90., 90. 
F,  0.0,  0.0,  0.0 
0 

0,  1,0,  1 
T,  4,  0.0,  0.0 
T,  5 
T 
7 
T 

T,  0.0 , 5.0,  V 
2,  '../eosdata' 

104,  104,  104,  20,  20,  1 


[title] 

[niter ,  var ,  varmin  ,  varmax] 

[restart ,  filnam] 

[xscale ,  tscalej 
[tmax  ,  maxcyl ,  safty] 

[lpexit,pexit,pratio.ncount,rnodestop,ntarget,ytmin,ytmax] 
[(slide  ,  ntarget ,  ytmin  ,  ytmax] 

[igeom  (0=Cartesian  ,  1  cylindrical)] 

[b1,b2,b3,b(0=r,1=t;  l,r,b,t)] 

[Ipburn  ,  ipgburn  ,  xdet ,  ydet] 

[Ifburn  ,  iffbum] 

['gauge] 

[Istres] 

[Ipfail] 

[iplot ,  tplot,  dtplot ;  splotdir] 

[ieos ,  eosfilnam] 

[material  EOS  numbers] 


Line  #1.  The  problem  title  is  input  on  the  first  line. 

Line  #2.  The  second  line  has  four  inputs  having  to  do  with  running  SMERF  in  an  iteration  mode.  A 
common  use  of  the  SMERF  code  is  to  determine  threshold  conditions  which  lead  to  the  detonation  of  an  acceptor 
explosive  material.  The  niter  iterations  will  be  done  to  converge  on  the  threshold  value  of  the  parameter  var  which 
separates  events  that  lead  to  a  detonation  from  those  that  do  not;  varmin  and  varmax  are  bounds  on  var.  If  niter 
is  zero,  then  no  iterations  will  be  performed  and  the  remaining  input  on  line  2  will  be  ignored  (it  should  still  be  in 
the  proper  format,  however).  The  parameter  var  is  a  two-character  string  (enclosed  in  single  quotes)  that  signifies 
what  parameter  will  be  iterated.  At  present,  there  are  four  possibilities.: 
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1.  var  =  ’pi':  The  shock  sensitivity  of  an  acceptor  explosive  material  will  be  iterated.  It  is  required  that 
the  Forest  Fire  bum  model  be  used  for  the  acceptor  material.  The  input  to  this  bum  model  is  the  straight 
line  on  the  Pop  plot  as  discussed  in  connection  with  equation  D4.1.  The  quantity  that  will  be  iterated  is  PI, 
the  pressure  (in  megabits  here  on  line  #2)  required  to  give  a  one-centimeter  run  distance.  The  material 
number  of  the  acceptor  material  is  given  by  the  parameter  iffburn  which  is  read  in  line  #11.  Further  input 
is  needed  to  determine  criteria  for  acceptor  explosive  detonation  or  no  detonation.  This  input  is  obtained 
from  line  #6.  If  the  problem  stop  time,  tmax,  or  stop  cycle,  maxcyl,  in  line  #  5  are  achieved,  then  the 
iteration  terminates,  so  these  values  have  to  be  set  sufficiently  high.  The  same  inputs  from  line  #5,  line  #6, 
and  line  #11  are  required  for  all  of  the  var  options. 

2.  var  =  'sc*.  The  scale  factors  that  are  normally  given  in  line  #4  will  be  iterated.  With  this  option,  the 
scale  factors  that  are  given  in  line  #4  will  be  read  by  SMERF  but  then  ignored.  This  option  automates  the 
computation  of  the  critical  diameter  of  an  explosive,  for  example.  Another  example  of  its  use  is  to 
determine  the  threshold  size  of  a  fragment  impacting  on  an  explosive.  It  is  assumed  in  the  iteration  that  a 
large  value  for  the  scale  factor  will  lead  to  a  detonation,  while  a  small  value  will  not. 

3.  var  =  'vy'.  The  initial  vertical  velocity  of  selected  materials  will  be  iterated.  The  initial  material 
geometry  is  given  in  the  matcor  file  described  later.  Included  in  the  geometry  description  is  the  initial 
velocity  of  each  material  region.  When  var  =  'vy',  each  non-zero  initial  vertical  velocity  is  set  to  the  same 
current  value  of  the  iterated  parameter.  This  option  is  used  to  find  the  threshold  velocity  for  detonation  for 
fragment  or  flyer  impact  on  explosive  targets. 

4.  var  =  Vx'.  The  initial  horizontal  (or  radial)  velocity  of  selected  materials  will  be  iterated.  This  option  is 
the  same  as  the  vertical  velocity  option  above. 

Line  #3.  If  there  is  to  be  a  restart  from  a  previous  run,  the  input  from  line  3  will  assign  the  value  true  to 
restart  and  the  name  of  the  restart  file  to  filnam.  Restart  files  are  dumps  of  the  COMMON  block  made  in  previous 
calculations.  The  frequency  of  dumps  during  a  run  is  controlled  by  the  outdata  file  described  in  Section  D,  Output 
Control  When  the  restart  option  is  set,  two  other  files  (in  addition  to  the  six  basic  input  files)  containing  history 
data  are  needed.  The  first  is  prestart,  which  contains  the  maximum  pressure  of  each  material  as  a  function  of  time. 
It  can  be  obtained  by  copying  pmatout  which  is  generated  by  each  SMERF  run.  Second,  if  the  lgauge  option  is 
true  (line  9  below),  then  the  file  grestart  is  required.  It  is  obtained  by  copying  the  file  gaugeout  that  results  from 
each  SMERF  run  whenever  lgauge=true.  It  should  be  noted  that  the  material  properties  can  be  changed  prior  to  a 
restart,  provided  that  the  material  property  that  is  changed  does  not  affect  the  mass  or  energy  of  the  material.  This  is 
done  by  editing  the  equation  of  state  data  file  (eosfilname  as  input  in  line  17).  It  should  also  be  noted  that  the 
restart  option  can  be  used  when  the  SMERF  code  is  in  the  iteration  mode  (line  2  above)  when  var  =  'pi',  and  the 
restart  file  was  created  prior  to  involvement  of  the  acceptor  explosive  in  the  event. 

Line  #4.  Scale  factors  that  multiply  all  input  spatial  dimensions  (xscale)  and  all  input  time  quantities 
(tscale)  are  set  in  this  line.  The  scale  factors  can  greatly  simplify  the  input  for  SMERF.  For  example,  the  length 
units  for  SMERF  are  centimeters,  but  it  is  sometimes  convenient  to  provide  input  in  inches  and  then  set  the  xscale 
to  2.54.  Sometimes  one  may  want  to  run  a  sequence  of  problems  where  the  only  change  is  the  size  of  the  objects 
involved.  It  is  easy  to  convert  the  input  for  one  size  problem  to  another  by  changing  the  xscale  factor.  In  this  case, 
one  should  also  change  the  time  scale  factor,  tscale,  by  the  same  amount. 

Line  #5.  The  problem  stop  time  (tmax),  stop  cycle  (maxcyl),  and  maximum  Courant  number  (safty)  are 
input  in  line  5.  The  Courant  number  governs  the  size  of  the  integration  time  step.  If  the  number  is  too  large,  then 
the  calculation  can  become  unstable,  and  it  will  terminate.  If  the  Courant  number  is  too  small,  then  the  computer 
may  become  excessive.  The  value  of  safty  can  range  from  about  0.4  to  0.5.  Sometimes  difficulties  will  be  obtained 
with  some  aspect  of  the  calculation,  and  the  program  will  terminate.  A  solution  might  be  to  repeat  the  calculation 
with  a  smaller  value  of  safty. 
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Line  #6.  Information  in  line  6  allows  execution  to  terminate  early,  depending  on  the  reaction  of  the 
explosive,  lpexit  =  true  stops  computation  if  detonation  is  detected  or  if  it  is  determined  that  no  detonation  will 
occur.  Detonation  is  assumed  if  the  pressure  in  the  explosive,  as  designated  by  the  material  number  iffburn 
described  below,  exceeds  pexit.  The  program  keeps  track  of  the  maximum  pressure  in  the  explosive.  Detonation 
failure  is  assumed  when  the  explosive  pressure  falls  by  a  factor  pratio.  When  detonation  or  failure  is  detected,  the 
program  terminates  ncount  cycles  later.  A  value  of  ncount  greater  than  zero  is  sometimes  useful  to  get  a  final  plot 
showing  development  of  the  detonation.  If  lpexit  =  true,  then  SMERF  will  terminate  in  one  of  two  modes.  For 
modestop=0,  the  calculation  stops  if  detonation  is  detected.  If  modestop^l,  the  run  stops  if  the  detonation  is 
produced  and  then  fails  to  propagate  in  the  x  direction  to  xstop  and  in  the  y  direction  to  ystop.  If  xstop  lies  outside 
of  the  mesh,  then  it  is  ignored.  However,  if  ystop  initially  lies  above  the  mesh,  the  mesh  may  eventually  be  shifted 
upward  (see  line  7)  to  include  ystop.  At  this  point,  ystop  becomes  effective. 


Line  #7.  If  lslide  =  true  in  this  line,  then  the  grid  will  be  moved  upward  to  follow  material  motion  such  as 
shock  waves  or  moving  projectiles.  If  motion  is  detected  at  the  top  of  the  grid,  then  another  row  of  elements  will  be 
added  to  the  top,  and  the  bottom  row  will  be  deleted.  This  can  save  enormous  amounts  of  computer  memory  and 
time  for  calculation  of  projectile  impact  onto  semi-infinite  explosives,  explosive  critical  diameter  experiments,  or 
self-forging  fragment  simulation.  Often  it  is  desired  to  simulate  impacts  with  targets  that  lie  outside  the  mesh.  This 
sliding  option  can  be  used  to  incorporate  a  target  with  material  number  ntarget  lying  between  vertical  coordinates 
ytmin  and  ytmax.  Whenever  this  option  is  used,  the  mesh  (as  specified  in  the  mshgen  data  file)  must  have  uniform 
spacing  in  the  vertical  direction. 

Line  #8.  The  coordinate  system  is  specified  in  line  8  as  Cartesian  (igeom-0)  or  cylindrical  (igeom=l). 

Line  #9.  Boundary  conditions  along  the  left,  right,  bottom,  and  top  (bl,b2,b3,b4)  of  the  grid  are  set  in  line 
9  and  can  be  either  reflective  (0)  or  transmissive  (1).  The  transmissive  boundary  allows  disturbances  and  material  to 
pass  out  of  the  computational  domain,  while  the  reflective  boundary  behaves  as  a  perfectly  reflecting  plane. 

Line  #10.  If  a  programmed  bum  is  desired,  then  lpburn=true  in  line  10.  The  remaining  three  parameters 
in  line  10  indicate  the  material  number  of  the  material  to  be  program  burned  (ipgburn)  and  the  point  of  initiation  of 
the  bum  (xdet,ydet).  The  material  chosen  to  bum  must  include  both  a  solid  and  gas  equation  of  state. 

Line  #11.  If  lfburn=true  in  line  11,  then  the  Forest  Fire  bum  is  enabled  for  any  material  assigned  this 
bum  option.  The  Forest  Fire  bum  is  useful  when  one  wishes  to  investigate  the  response  of  an  explosive  to  a  shock 
stimulus.  The  second  entry  in  line  1 1,  iffburn,  is  the  material  number  of  the  “acceptor”  explosive.  This  number  is 
required  when  Ipexit=true  in  line  6  above.  In  addition,  when  Ifburn=true,  the  maximum  pressure  from  material 
iffburn  is  included  in  the  printed  output  from  the  calculation.  Otherwise,  the  maximum  pressure  of  all  materials  is 
printed. 


Line  #12.  If  Igauge=true  in  line  12,  then  the  time,  pressure,  and  speed  at  locations  specified  in  input  file 
gaugein  are  written  each  cycle  to  output  file  gaugeout. 

Line  #13.  If  material  strength  is  to  be  included  in  the  problem,  Istres=true  in  line  13. 

Line  #14.  Tensile  failure  of  materials  is  allowed  if  Ipfail=true  in  line  14. 

Line  #15.  Line  15  contains  switches  for  plotting  data.  When  Iplot=true  and  the  simulation  time  exceeds 
tplot,  the  data  for  pressure  contour  plots  is  dumped  every  dtplot  microseconds  into  directory  splotdir. 
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Line  #16.  The  integer  ieos  is  a  variable  that  is  read  but  not  used  by  the  SMERF  program.  The  file  name 
for  the  appropriate  equation-of-state  data  package  is  eosfil. 

Line  #17.  The  last  line  of  setup  contains  the  equation-of-state  numbers  of  the  materials  in  the  problem. 
These  numbers  identify  materials  found  in  the  equation-of-state  data  file  eosfilnam. 


B.  THE  MESH  (The  mshgen  File) 

The  computational  mesh  is  generated  using  data  in  input  file  mshgen  which  supplies  the  coordinates  for 
the  entire  mesh  (xj,X2,yj,y2)  and  a  finely  zoned  subgrid  (xsg|,xsg2,ysgj,ysg2)  with  constant  cell  size  (dxs,dys). 
Zones  grow  outwardly  from  the  subgrid  by  the  factors  xfac,yfac.  The  number  of  materials  in  the  problem  (nmat)  is 
input  in  the  first  line  of  mshgen.  The  example  problem  discussed  in  Section  IV  uses  a  simple  mesh  with  constant 
0.2-cm  square  cells  throughout.  The  mshgen  file  for  it  looks  like  the  one  shown  in  Table  3. 


TABLE  3.  mshaen  Input  File  for  the  Bomb 
Stack  Example  Problem. 

6  [nmat] 

0.00,  45.00  [xl  .  x2] 

0.00,  45.00  [xsgl  .  xsg2] 

0.20,  1.00  [dxs  ,  xfac] 

0.00,  45.00  [yl  ,  y2] 

0  00,45.00  [ysgl  ,  ysg2] 

0.20,  1.00  [dys  ,  yfac] 


C.  MATERIAL  ZONING  (The  matcor  File) 

The  last  line  in  setup  contains  the  equation-of-state  numbers  of  the  different  materials  in  the  problem. 
These  numbers  are  used  to  assign  equation-of-state  data  to  the  materials.  The  location  of  each  material  in  the  mesh 
is  determined  by  the  data  supplied  in  the  file  matcor.  This  geometry  data  is  in  the  form  of  a  sequence  of  ordered 
points  and  arcs.  The  data  read  in  matcor  specifies  the  boundary  of  a  material  body.  The  data  is  in  the  form  of  an 
integer  (itype)  followed  by  a  coordinate  pair  defining  a  point.  The  options  available  for  itype  are  listed  in  Table  4 
and  illustrated  in  Figure  5.  For  some  values  of  itype,  the  coordinate  represents  a  vertex  on  a  polygon.  For  other 
values  of  itype,  the  coordinates  represent  the  center  of  an  arc.  In  this  case,  the  counterclockwise  arc  length  in 
degrees  is  then  read  by  SMERF  in  the  following  line.  The  first  entry  must  have  itype=l,  and  the  coordinate  pan- 
then  represents  a  starting  point  on  the  material  boundary.  The  sequence  is  terminated  by  the  coordinate  pair  (nexit, 
0.0),  where  nexit  is  a  number  greater  than  1.0e+23. 
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The  next  line  in  matcor  contains  the  initial  state  for  the  material  object  that  was  just  defmed.  The  six 
entries  on  this  line  give  the  material  number,  starting  pressure  (megabars),  temperature  (degrees  Kelvin),  horizontal 
velocity  (cm/ps),  vertical  velocity  (cm/ps),  and  history  variable.  If  the  material  has  a  gas  equation  of  state,  then  the 
history  variable  is  the  mass  fraction  of  the  explosive  gas  products.  If  the  material  is  porous,  then  the  history  variable 
is  the  initial  porous  volume  divided  by  the  specific  volume  of  the  solid  matrix. 

Sometimes  a  material  object  will  be  used  repeatedly  at  different  places  on  the  mesh  or  it  may  be  more 
convenient  to  describe  it  at  one  location  and  then  rotate  and  move  it  to  a  different  location.  One  may  define  a 
material  geometry  prototype  by  entering  itype=0  followed  by  a  coordinate  pair.  The  coordinate  pair  in  this  case 
represents  a  reference  point  on  the  prototype  geometry.  The  material  geometry  boundary  is  then  entered  in  the  same 
manner  as  explained  above  until  terminated  by  the  coordinate  pair  (nexit,  0.0).  The  prototype  geometry  data  is 
saved  in  a  buffer.  At  any  point  during  the  matcor  data  input,  the  prototype  geometry  can  be  used  to  construct  a 
material  object  on  the  mesh.  This  is  done  by  entering  a  line  with  itype=-l  and  a  destination  coordinate  pair 
followed  by  an  angle  on  the  next  line.  The  prototype  geometry  region  is  rotated  in  the  counterclockwise  direction 
about  the  reference  point,  and  then  the  entire  object  is  translated  until  the  reference  point  coincides  with  the 
destination  point.  A  material  properties  line  is  then  read  in  the  usual  format.  The  prototype  geometry  can  be  used 
as  many  times  as  desired  to  form  many  material  objects.  Each  material  object  can  have  different  materials  and 
different  initial  conditions.  Only  one  prototype  can  exist  at  a  time.  If  data  for  a  second  prototype  geometry  is 
entered,  it  replaces  the  previous  one. 


TABLE  4.  Options  for  Connecting  the  Points  Specified  in  File  matcor. 


itype 

point 

relative  to 

last  point 

r 

(x,y) 

(0,0) 

line  ; 

■  2 

(r,0) 

(0,0) 

line 

I-  3.. ; 

(*cM 

(0,0) 

arc 

ti 

lilt  §0Wm 

last  pt. 

line  : 

12 

m 

last  pt. 

line 

13 

(rC’0c) 

last  pt. 

arc 
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Figure  5.  Different  Ways  to  Connect  Two  Points  in  the  SMERF  Setup.  Material  regions  are 
defined  on  the  mesh  by  using  one  or  more  of  these  options  to  generate  geometrical  shapes. 
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The  matcor  file  used  in  the  bomb  stack  example  problem  in  Section  IV  is  shown  in  Table  5.  This  example 
uses  the  prototype  geometry  option.  A  prototype  bomb  case  is  defmed  at  the  coordinate  origin,  for  example,  and 
then  it  is  moved  to  construct  the  donor  and  three  acceptor  bomb  cases. 

TABLE  5.  matcor  Input  File  for  the  Bomb  Stack  Example  Problem. 


1,  -10000. ,  -1G0GO. 

1,  -10000.,  +1OG0G. 

1,  +10000., -10000. 

1,  +10000., -10000. 

1,  S.e+23, 0. 

6, I.Qe-6,  300.,  0.0,  0.0,  0.0 
-1,0.0,  0.0 
1,0.0,  13.65 
3,  0.0,  0.0 
360. 

1,  5.0e+23,  0. 

0,0.0,  0.0  0.0 

1, 1.0e-6,  300.,  0.0,  0.0,  0.0 
0,  28.57,  0.0 
0.0 

2,  T.Oe-6,  300.  0.0,0.0,  0.0 
0,  0.0,  28.57 

0.0  :  " 

2,  1 .0e-6,  300.,  0.0 ,  0.0,  0.0 
0,  28.57,  28.57 

0.0 

3,  T.Qe-6, 300.,  0.0,  0.0,  0.0 
-1,0.0,00 
1,0.0,12.65 

3,  0.0 , 0.0 
360. 

1 , 5.0e+23,  0. 

0,  0.0 ,0.0 
0.0 

4,  1.0e-6,  300.,  0.0 ,0.0, 0.0 
0,  28.57 , 0.0 

0:0 

5, 1.0e-6, 300.  ,0.0,0.0,00 
0,  0.0 , 28.57 
0.0 

5, 1  Oe-6,  300. ,  0.0  ,0.0  ,  0.0 
0,  28.57 , 28.57 
OO 

5, 1.0e-6 , 300.  ,  0.0 ,0.0,00 


Air 

Case  Prototype 
Case  1,1 

:  Case  1,2 

Case  2,1 

Case  2,2 

HE  Prototype 
HE  1,1 

HE  1,2 

HE  2,1 

HE  22 


There  are  several  important  rules  concerning  input  of  materials  into  the  mesh: 


Rule  #1.  If  material  regions  overlap,  the  last  region  specified  is  used.  For  example,  a  rectangular 
steel  plate  with  a  circular  hole  would  be  constructed  by  first  specifying  the  steel  rectangle 
followed  by  a  circle  of  air. 
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Rule  #2.  Every  computational  zone  must  contain  a  material  Air  is  often  the  best  material  to  use 
to  fill  the  mesh  and  play  the  part  of  vacuum. 

Rule  #3.  Air  must  be  the  last  material  number  as  input  in  data  file  setup.  Air  must  be  included  in 
all  problems.  It  must  be  the  last  material. 

Rule  #4.  Material  regions  can  extend  beyond  the  mesh  boundary.  SMERF  truncates  the  regions 
to  fit  within  the  mesh. 

Rule  #5.  Material  regions  (polygons)  are  automatically  closed.  For  example,  a  triangle  can  be 
specified  by  three  entries. 


D.  OUTPUT  CONTROL  (The  outdata  and  gaugein  Files) 

Three  types  of  output  are  available  from  SMERF:  (1)  binary  COMMON  dumps,  which  can  be  used  for 
restarts  or  as  input  to  graphics  packages;  (2)  material  maps  written  to  the  screen;  (3)  contour  plots;  and  (4) 
Lagrangian  time  histories  for  pressure  and  particle  speed.  The  timing  of  the  COMMON  dumps,  material  maps,  and 
contour  data  is  set  in  outdata  where  a  ‘D’  or  ‘d’  (dump  COMMON),  an  ‘M’  or  ‘m’  (material  map),  or  a  ‘P’  or  ‘p’ 
(contour  plot)  followed  by  the  output  time  is  specified.  The  material  map  is  a  nx  x  ny  matrix  of  symbols  indicating 
the  location  of  the  various  materials  in  the  computational  mesh  at  a  given  time.  This  map  is  a  simple  way  to  view 
the  state  of  the  calculation  as  it  is  running.  It  is  written  to  standard  output.  For  the  example  problem  discussed  in 
Section  IV,  the  file  outdata  looks  like  the  example  in  Table  6. 

TABLE  6.  outdata  File  for  the  Bomb 
Stack  Example  Problem. 


W ,  0.0 

M,  1.5 
•D* ,  2.0 

•M1 ,  2.5 
F,  3.0 


The  Lagrangian  history  of  the  flow  is  often  required.  The  initial  locations  of  the  tracer  particles  are  specified  in 
gaugein.  The  first  entry  is  the  number  of  gauge  points,  with  coordinates  on  the  lines  following.  A  simple  example 
file  containing  two  points  is  given  below. 


TABLE  7.  A  gaugein  File  With  Two 
Lagrangian  History  Points. 


2 

0.000 ,5.000 
2.335 ,8.990 
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E.  EQUATION  OF  STATE  CONSTANTS  (The  eosdata  File) 

Equation-of-state  constants  for  materials  are  found  in  an  equation-of-state  data  file.  The  library  supplied 
with  the  code  is  called  eosdata.  The  first  line  for  each  material  contains  two  entries,  the  equation-of-state  number 
assigned  to  that  material  and  the  name  of  the  material.  The  second  line  contains  five  entries  which  assign  specific 
models  to  the  material:  nmat  (material  strength),  neoss  (solid/reactant  eos),  neosg  (gas/product  eos),  nburn  (high 
explosive  bum),  npor  (porosity).  The  solid/reactant  eos  category  designates  the  reactant  or  original  state  of  the 
material  (whether  it  is  in  solid  or  gas  form),  thus  the  inclusion  of  the  “perfect  gas”  eos  in  this  category.  The 
gas/product  eos  category  is  reserved  for  the  equations  of  state  of  the  products  of  chemical  reaction  in  explosive 
materials. 

The  options  for  the  constitutive  model  categories  are  described  in  Table  8  (in  reality,  a  group  of  four 
tables).  Inspection  of  the  table  shows  that  the  units  of  the  input  parameters  are  not  entirely  consistent.  In  the  early 
development  of  the  SMERF  equation-of-state  subroutines,  material  data  was  obtained  from  a  variety  of  databases 
which  sometimes  used  different  units.  In  particular,  the  databases  used  were  for  the  DYNA2D  code  from  the 
Lawrence  Livermore  National  Laboratory  and  the  2DE  code  from  the  Los  Alamos  National  Laboratory.  Not  all  of 
the  features  from  these  codes  were  implemented  so  that  some  of  the  constants  in  the  present  database  are  not  used. 
The  entries  in  eosdata  for  the  materials  used  in  the  example  problems  are  listed  in  Table  9. 


TABLE  8.  Definition  of  the  Equation  of  State  Constants  Read  in  the  eosdata  File. 


NMAT 

MO  DEL/E QNS 

NAME 

UNITS/COMMENTS 

0 

None 

Only  for  perfect  gas 

9 

Fluid 

Po 

initial  Density 

g/cm'3 

Pc 

Spall  Pressure 

Mb  (must  be  negative) 

10 

Elastic- 

Po 

Initial  Density 

g/cm'3 

Perfectly  Plastic 

ft 

Shear  Modulus 

Mb  1 

Y 

Yield  Strength 

Mb 

E 

Hardening  Modulus 

not- used 

Pc 

Spall  Pressure 

Mb  (must  be  negative) 

:  ■  -  A,  ;;  :V': 

Yield  Constant 

not  used 

a2 

Yield  Constant 

not  used 

Ps 

Yield  Constant 

not  used 

NEOSS 

MODEL/EQNS 

SYMBOL 

NAME 

UNfTS/COMMENTS 

2 

Perfect  Gas 

Y 

Ratio  of  Sp.  Heats 

dimensionless 

Eq.  Cl  .2 

Q 

Energy 

Mb 

4 

Gruneisen 

Hugoniot 

cm/us 

Eq.  Cl  3 

S1,S2,S3 

Coefficients 

.  WMi 

Gruneisen  Coef. 

sa 

Vol  Dependence  of 

dimensionless 

120 

HOM 

A,B,C,D,E 

Coefs.  in  Eq.  Cl. 5c 

Eq.  C  E5 

K,L,MfN,0 

Coefs.  in  Eq.  C1.5e 

Q.R,S,T,U 

Coefs.  in  Eq.  G1.5d 

Cv 

i  Specific  Heat 

cal/g 

Z 

Energy  Offset 

Mbar-  cnrVg 
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TABLE  8  (Cont’d).  Definition  of  the  Equation  of  State  Constants  Read  in  the  eosdata  File. 


NMAT 

MODEL/EQNS 

SYMBOL 

NAME 

UNITS/COMMENTS 

3 

Arrhenius 

E 

Activation  Energy 

kcal 

z 

Frequency  Factor 

sec'1 

7 

ignition  &  Growth 

Eqs.  D3.2-D3.4 

Sffti  ill  ekflMfe 

ii 

Forest  Fire 

n 

not  used 

Eq.  D4.1 

tiiii 

Pop-Plot  Point 

Kb  -  note  the  units! 

-s 

Pop-Plot  Slope 

must  be  negative 

p2 

not  used 

6 

JWL 

A ,  B 

Coefficients 

Mb 

Eq.  Cl  .6 

Rl  ,  ^2 

Exponents 

dimensionless 

wCv 

Gruneisen  x  Sp  Heat 

Mb-  cnf/g-K 

Cv 

Specific  Heat 

Mb-  cm7g-K 

E0  . 

Energy 

Mb-  cm% 

100 

HOM 

Ci-S, 

C  is  cm/ms 

Eq.  Cl. 4 

not  used 

Gruneisen  Coef. 

dimensionless 

m  :  ■ 

Specific  Heat 

cal/g 

VO 

not  used 

a  ' 

Thermal  Expansion; 

crri/K 

PS.US,T0 

not  used 

NPOR 

MODEL/EQNS 

SYMBOL 

NAME 

UNITS/COMMENTS 

0 

None 

3 

ii  P-a 

ns 

•  *  :  Pn 

Initial  Density 

gf  cm3 

Elastic  Sound  Speed 

cm/ps 

(Vci,Pci,i=l.ns) 

Compaction  Curve 
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TABLE  9.  Entries  in  eosdata  for  Materials  Used  in  the  Bomb  Stack  Example  Problem. 


1, . ’AIR' . 

1,2,0,  0,0, 

1.4,  1.0e-5,  0.0, 

20,  'AFX-1100 ', 

10,110,6,11,0, 

1.492,  .0500,  .001,  0.0,  -0.01, 

0,  0,  1, 

.206,  2.16,  .01,  0.00,  0.00, 

1.50,  .259,  .670,  5.0e-5,  0.0, 

0.0,  300.0, 

3.91,  0.0252,  4,71,  1.18,  2:0e-6,  1e-5,  0.054, 

6,  84.4,  -2.216,  0.0,  [  PI  =  84.4kbar  ] 

104,  'HOM  STEEL', 

10,110,0,0,0, 

7.9,  0.769,  .010,  0.0,  -0.020, 

0,  0,  1, 

4.580E-01 ,  1.510E-00,  6.500E-02, 

0.000E-00,  0.000E-00, 

2.020E-00,  1.070E-01,  1.266E-01, 

1.170E-05,1  290E-01, 

1.200E-01,  3.000E+02, 

117,  'HOMCOMP-B', 

10,  110,  120,  11, 0, 

1.715,  .0500,  .002,  0.0,  -0.01, 

0,  0,  1, 

.280,  1  63,  .001,  0.00,  0.0, 

1.5,  .259,  .583,  5.0E-5,  0.0, 

0.0,  300, 

-35029893E+0 1 ,  -.24367332 E+01,  .267361 84E+00, 
-.24236604E-01 ,  .66789557E-03, 

-.  12981 939E+01,  .31351829E+00,  46863277E-01 , 
.32738625E-02,  .85793089E-04, 

. 75364753 E+01,  -.483730 16E+00,  10225829E+00, 

-.1529242 IE-01,  .76422916E-03, 

5,  .17530137, 

6.53.9,  -1.501  ,0.0,  { P1=53.89  ] 
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F.  BASIC  OUTPUTS 

There  is  no  post-processor  for  the  SMERF  code.  Certain  programs  have  been  written  by  the  authors  of 
SMERF  for  graphical  presentation  of  the  output,  but  the  programs  are  machine-dependent  and,  thus,  are  not 
normally  distributed.  These  programs  read  data  in  output  files  produced  by  SMERF.  In  this  section,  the  contents 
and  format  of  these  files  are  explained.  The  names  of  subroutines  in  SMERF  where  these  output  files  are  produced 
are  given  so  that  the  interested  code  user  can  modify  them. 

COMMON  dumps  are  both  written  and  read  in  subroutine  RW.  This  output  is  unformatted  and  contains 
essentially  all  the  information  regarding  the  state  of  the  calculation  at  a  particular  time  step.  These  files  can  be 
extremely  large  and  are  used  to  restart  SMERF  calculations  as  explained  in  the  setup  file  description. 

Because  of  the  large  size  of  the  COMMON  dumps,  another  means  of  producing  data  for  two-dimensional 
contour  plots  is  given.  These  output  files  are  called  splots.  In  the  present  version  of  the  SMERF  code,  sufficient 
data  is  output  for  the  production  of  contour  plots  of  pressure  and  the  history  variable  (bum  fraction  or  porosity), 
together  with  material  boundaries.  The  timing  of  the  output  of  these  files  can  be  controlled  through  the  outdata  file 
using  the  "P"  option  as  described  above,  splots  can  also  be  generated  at  regular  intervals  by  means  of  input  data  in 
the  setup,  also  described  above.  The  directoiy  where  the  splots  are  written  is  also  specified  in  the  setup.  Names  of 
these  files  are  assigned  according  to  the  formula  splotxxxx  where  xxxx  is  the  sequence  number  of  the  files  starting 
at  0001.  They  are  formatted  and  can,  therefore,  be  viewed  and  edited.  Another  reason  for  formatting  the  files  is 
that  the  primary  SMERF  calculations  are  sometimes  performed  on  a  computer  which  is  different  from  the  one  used 
for  post-processing.  Formatted  files  can  be  transferred  over  the  network  without  any  word  length  conversions.  The 
writing  of  the  files  is  done  in  subroutine  splot.  These  files  are  most  often  used  to  produce  animated  sequences  of 
contour  plots  of  the  simulation. 

If  the  gauge  option  was  enabled  in  the  setup,  then  an  output  file  gaugeout  is  produced.  It  contains 
formatted  data  representing  the  time  histories  of  pressure  and  particle  speed  at  specified  gauge  points  which  move 
with  the  material.  This  file  is  produced  in  subroutine  GAUGE. 


The  final  output  file  is  pmatout.  It  contains  time  histories  of  the  maximum  pressure  for  each  material  in  the 
SMERF  calculation.  This  formatted  file  is  produced  in  subroutine  ASTOP. 
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IV.  EXAMPLE  PROBLEMS 


A.  FRAGMENT  IMPACT 


The  first  example  problem  concerns  shock  initiation  of  an  explosive  from  fragment  impact.  This  is  "demo” 
problem  #1  supplied  with  the  code  on  installation.  The  "fragment"  is  taken  to  be  a  steel  pellet  with  a  diameter  of 
1.00  cm,  height  of  0.20  cm,  and  speed  of  0.12  cm/ps.  The  explosive  is  a  cylinder  of  bare  COMP-B  with  a  diameter 
of  40  cm  and  a  height  of  20  cm.  A  schematic  of  the  setup  is  shown  in  Figure  6. 


This  simulation  employed  a  uniform  rectangular  subgrid  (0.50  cm  x  2.00  cm)  embedded  within  a  1.00  cm 
x  2.00  cm  computational  domain.  The  grid  cells  comprising  the  fine  scale  subgrid  were  0.03  cm  square.  These 
zones  were  allowed  to  expand  by  3%  in  the  radial  direction  for  each  column  outside  the  subgrid,  but  the  vertical  cell 
dimension  was  not  allowed  to  grow.  In  this  way,  30  x  66  total  grid  cells  were  generated  for  use  in  the  simulation. 
Input  files  for  this  example  problem  are  shown  in  Tables  10,  11,  and  12.  Results  of  the  calculation  are  shown  in 
Figure  7. 


f 


j 

1.5  cm 


i 


V 


t 

0.2  cm 

I 


(1.00,1.5o) 


Comp-B 


Steel 

6.50  cm 


0.12  cm/us 


Air 


0.00, -0. 


50) 


FIGURE  6.  Schematic  of  the  Setup  for 
Shock  Initiation  by  Fragment  Impact. 
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TABLE  10.  setup  Input  File  for  the  Fragment  Impact  Example  Problem. 


'Fragment  Impact' 

0,  ’uy‘,  0.1,  0.14 
F.'restartfile' 

1.0, 1.0 

3.0,  500,  0.4 

T,  .2,  .3,  100, 0,90., 90. 

T,  0.0,  20.0,  20.5 
1 

0,1,  iii 

F,  0,  0.0,  0.0 
T,  1 

m. 

■T  '■  - 

;:T '  ’ 

T,  0.0 , 0.50, 

2,  '../eosdata' 
117,104,1 


[title] 

[niter ,  var ,  varmin  ,  varmax] 
[restrt ,  filnam] 

[xscaie ,  tscale] 

[tmax  ,  maxcyl ,  safty] 

[Ipexit  ,  pexit ,  pratio,  ncount, 
modestop ,  xstop ,  ystop] 
[Islide  ,  ntarget ,  ytmin  ,  ytmax] 
[igeom] 

[b1,b2,b3,b] 

[Ipburn  ,  ipgburn  ,  xdet ,  ydet] 
[Ifburn  ,  iffbum] 

[igauge] 

[Istres] 

[Ipfail] 

[(plot ,  tplot,  dtplot ,  splotdir] 
[ieos ,  eosfil] 

[material  EOS  numbers] 
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TABLE  11.  mshaen  Input  File  for  the  Fragment  Impact  Example  Problem. 


0.00, 1.00 
0.00 , 0.50 
0.03,  1.03 
-0.05 , 1.50 
-0.05 ,1.50 
0.03 , 1.00 


[nmat] 

[X^] 

[xsg  l  ,  xsg2J 
[dxs ,  xfac] 
M.YJ 
[ysgl ,  ysg2] 
[dys ,  yfac] 


TABLE  12.  matcor  Input  File  for  the  Fragment  Impact  Example  Problem. 


IV 

1, 

1, 

1, 

1, 


0.00  ,  -60.00 
25.0  ,  -60.00 
25.0,  60.00 
0.00,  60.00 
5.e+23 , 0. 


3,  1.0e-6.  300. ,  0.0,  0.0,  0.0  AIR 


1, 
1, 
IV 
1 , 
1, 


0.00,  0,00 
20.0,  0.00 
20.0  ,  20.00 
0.00  ,  20.00 
5.e+23 , 0. 


1,  1.Ge-6,  300.,  0.0,  0.0,  0.0  COMP-B 


IV 

1, 

iv 

1, 

1, 


0.00  ,  -0.20 
0.50  ,  -0.20 
0.50,  0.00 
0.00  ,  0.00 
5,6+23 ,0. 


2, 1  .Oe-6,  300.,  0.0,  0.12,  0.0  STEEL 
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Figure  7c.  Early  detonation  at  2.0  ms.  The  maximum  Figure  7d.  Mature  detonation  wave  in  the  Comp-B  at 

pressure  is  250  Kb.  2.5  ps. 


TheUpressure^^Econtolureda,inllg«^tscale?0<One,'set'ol  solid  ^acl^'lines^s^the^O^reactlo^extent 

contour  and  the  other  delineates  material  (air-fragment-explosive)  boundaries.  The  fragment 
was  modeled  as  a  steel  pellet  1.0  cm  in  diameter  and  0.2  cm  thick  with  speed  of  1.2  km/s. 
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B.  SYMPATHETIC  DETONATION  IN  BOMB  STACKS 


The  second  example  problem  is  sympathetic  detonation  in  a  stack  of  four  bombs.  This  is  "demo"  problem 
#4  supplied  with  the  code  on  installation.  The  bombs  are  modeled  as  infinitely  long,  13.65-cm-radius  cylinders  with 
1.00-cm-thick  steel  walls.  The  explosive  fill  is  AFX-1100,  an  experimental  Air  Force  explosive  that  has  been 
extensively  characterized.  Detonation  of  the  donor  bomb  is  modeled  using  the  programmed  bum  option  with 
detonation  initiated  at  the  center  of  the  bomb.  The  Forest  Fire  bum  option  is  used  to  simulate  shock-induced 
reaction  in  the  acceptor  bombs.  A  schematic  diagram  of  the  setup  is  shown  in  Figure  8.  The  input  files  for  this 
problem  are  not  tabulated  here.  They  were  used  as  examples  in  Section  III,  which  discussed  code  input,  and  can  be 
found  there. 

Sympathetic  detonation  in  a  simple  two-body  configuration  of  these  bombs  has  been  tested  [15].  The 
bombs  are  placed  side  by  side  and  the  donor  bomb  is  detonated.  With  the  AFX- 1 100  fill,  the  acceptor  bomb  was 
never  observed  to  detonate,  regardless  of  the  separation  distance  between  the  two  bombs.  The  bombs  were  also 
tested  in  a  stack  configuration.  Four  bombs  were  arranged  in  a  rectangular  stack  with  a  1.27-cm  gap  between 
adjoining  bombs.  It  was  found  experimentally  that  sympathetic  detonation  occurred  in  the  acceptor  bomb  located 
diagonally  from  the  donor  bomb.  Figures  9a  through  9d  show  pressure  contours  in  the  simulated  bomb  stack  at  4 
different  times.  The  first  figure  shows  programmed  detonation  in  the  donor  bomb.  The  second  figure  shows  the 
donor  case  accelerating  toward  the  diagonal  weapon.  The  third  figure  shows  the  onset  of  burning  in  the  diagonal 
acceptor.  Sympathetic  detonation  is  clearly  seen  in  the  last  figure. 

In  this  calculation  it  is  important  to  note  that  transition  to  detonation  is  a  natural  result  of  the  combustion  of 
the  explosive  which  is  modeled  solely  on  the  basis  of  independent  small-scale  tests.  There  are  no  additional 
switches  that  turn  detonation  off  or  on.  The  result  is  a  true  prediction. 


FIGURE  8.  Schematic  of  the  Setup  for  Simulation  of 
Sympathetic  Detonation  in  the  Mk-82  Bomb  Stack 
With  AFX-1100  Explosive  Fill. 
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Figure  9a.  Programmed  bum  in  the  donor 
bomb  at  15  ps. 


Figure  9b.  Donor  case  accelerating 
towards  the  diagonal  acceptor  at  95  ps. 


Figure  9c.  Onset  of  bum  in  the  addeptor 
due  to  impact  shock  at  135  ps. 


Figure  9d.  Sympathetic  detonation  in  the 
diagonal  weapon  at  160  ps. 


FIGURE  9.  SMERF  Simulation  of  Sympathetic  Detonation  in  the  Mk-82  Bomb  Stack.  Steel  case  material  is 
identified  with  solid  black  lines  and  the  pressure  field  is  contoured  in  gray  scale.  A  ‘programmed  burn’  was 
used  to  detonate  the  donor  bomb  and  a  ‘Forest  Fire*  burn  was  used  to  compute  the  reaction  in  the  three 
acceptor  bombs.  The  calculation  predicts  detonation  only  in  the  diagonal  weapon  due  to  impact  from  the 
expanding  donor  case.  The  solid  contour  within  the  burn  region  in  figure  c  is  the  50%  reaction  extent 
demarkation.  The  explosive  fill  used  was  AFX-1100. 
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